Collapse of magnetized hypermassive neutron stars in general relativity 
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Hypermassive neutron stars (HMNSs) - equilibrium configurations supported against collapse by 
rapid differential rotation - are possible transient remnants of binary neutron star mergers. Using 
newly developed codes for magnetohydrodynamic simulations in dynamical spacetimes, we are able 
to track the evolution of a magnetized HMNS in full general relativity for the first time. We find that 
secular angular momentum transport due to magnetic braking and the magnetorotational instability 
results in the collapse of an HMNS to a rotating black hole, accompanied by a gravitational wave 
burst. The nascent black hole is surrounded by a hot, massive torus undergoing quasistationary 
accretion and a collimated magnetic field. This scenario suggests that HMNS collapse is a possible 
candidate for the central engine of short gamma-ray bursts. 
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Hypermassive neutron stars (HMNSs) figure promi- 
nently in several relativistic astrophysical systems of cur- 
rent interest. Mass limits for nonrotating stars [the 
Oppenheimer-Volkoff (OV) limit] and for rigidly rotat- 
ing stars (the supramassive limit, which is only about 
20% larger) can be significantly exceeded by the pres- 
ence of differential rotation [l| . Stars with masses greater 
than the supramassive limit are called hypermassive 
stars. Mergers of binary neutron stars (BNSs) can lead 
to the formation of HMNS remnants. This possibility 
was explored in Newtonian and post-Newtonian simula- 
tions 0, 12] , and in full general relativity Q| . The latest 
relativistic BNS merger simulations with realistic equa- 
tions of state (EOS) [| confirm that HMNS formation 
is indeed a possible outcome. HMNSs can also result 
from core collapse of rotating massive stars, since rapid 
differential rotation can develop during the collapse. 

Differentially rotating stars approach rigid rotation via 
transport of angular momentum on secular time scales. 
HMNSs, however, cannot settle down to rigidly rotating 
equilibria since their masses exceed the maximum allowed 
by uniform rotation. Thus, 'delayed' collapse to a black 
hole, and possibly mass loss, will follow transport of an- 
gular momentum from the inner to the outer regions. 
Previous calculations of HMNS collapse have focused on 
viscous angular momentum transport 6] and angular mo- 
mentum loss due to gravitational radiation p|. In this 
Letter, we demonstrate black hole formation induced by 
seed magnetic fields in HMNSs. 

In any highly conducting astrophysical plasma, a 
frozen-in magnetic field can be amplified appreciably 
by gas compression or shear (e.g. differential rotation). 
Even when an initial seed magnetic field is weak, the field 
can grow to influence significantly the system dynamics. 
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There are at least two distinct effects which amplify the 
magnetic field in a HMNS: magnetic winding and the 
magnetorotational instability (MRI) [Lsl IT^ | . Numerical 
simulations using a general relativistic magnetohydrody- 
namics (GRMHD) code are required to follow this growth 
and determine the consequences. The key subtlety is that 
the wavelength of the fastest-growing MRI mode must 
be well resolved on the computational grid. Since this 
wavelength is proportional to the magnetic field strength, 
it becomes very difficult to resolve for small seed fields. 
However, the simulations reported here succeed in resolv- 
ing MRI. 

New computational tools now make long-term numeri- 
cal evolutions of relativistic magnetized HMNSs possible 
for the first time. Two groups have independently devel- 
oped codes for evolving MHD fluids in strong-field, dy- 
namical spacetimes []J0| (see also @). These codes solve 
the Einstein-Maxwell-MHD system of coupled equations, 
both in axisymmetry and in 3+1 dimensions, essentially 
without approximation. Both codes evolve the space- 
time metric using the BSSN formulation [lOj and em- 
ploy conservative shock-capturing schemes to integrate 
the GRMHD equations. Multiple tests have been per- 
formed with these codes, including MHD shocks, MHD 
wave propagation, magnetized Bondi accretion, MHD 
waves induced by gravitational waves, and magnetized 
accretion onto a neutron star. 

To study the effects of magnetic fields on HMNSs, we 
first construct initial data assuming a T — 2 polytropic 
EOS, P = Kp r , where P, K, and p are the pressure, 
polytropic constant, and rest-mass density. Henceforth, 
we adopt units such that K = c = G = 1, where c is 
the speed of light and G is the gravitational constant. In 
these units, the maximum allowed baryon masses Mq of 
nonrotating and o f rig idly rotating stars are 0.180 and 
0.207, respectively [Tlj . 
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rotation. In this paper, A is set equal to the coordinate 
equatorial radius R, giving a value of ~ 1/3 for the ratio 
of equatorial to central fi. 

We model a HMNS with Mq — 0.303, maximum den- 
sity pmax = 0.0668, and angular momentum parameter 
J/M 2 = 1.00 [Hj. The Arnowitt-Deser-Misner (ADM) 
mass is given by M = 0.279, which is about 70% larger 
than the OV limit. Because of rapid differential rota- 
tion, the shape of this star is highly flattened (see the 
first panel of Fig. 2). 

During the evolution, we adopt a T-law (adiabatic) 
EOS P = (r - l)pe with V = 2. Here, e denotes 
the specific internal energy. Before evolving, we add a 
weak poloidal magnetic field to the equilibrium model 
by introducing a vector potential of the form A v — 
n7 2 maxL4{,(P — P cut ),0], where P cut is 4% of the maxi- 
mum pressure, Aj, is a constant that determines the initial 
strength of the magnetic field, and w is the cylindrical 
radius. (A similar form of vector potential has been used 
in other MHD simulations 14].) We characterize the 
strength of the initial magnetic field by C = max(6 2 /P), 
the maximum value on the grid of the ratio of the mag- 
netic energy density b 2 to the pressure. Several val- 
ues of At are chosen to yield the following values of C: 
1.26 x 10~ 3 , 2.47 x 10~ 3 ', 4.88 x 10~ 3 , and 9.80 x 10~ 3 . 
We have verified that these small initial magnetic fields 
introduce negligible violations of the Hamiltonian and 
momentum constraints. Note that C oc v\, where va is 
the Alfven speed. Thus, C~ x / 2 cx uT 1 is proportional 
to the Alfven time tA — R/va- If the evolution time 
scale is determined by the Alfven time, a scaling rela- 
tion in the evolution should hold. This relation may not 
hold for all C due to MRI, which grows exponentially at 
a rate independent of tA- Comparing simulations with 
the values of C quoted above, we indeed find that, if the 
time is rescaled as C~ 1 ^ 2 t, the results for different values 
of C are approximately the same. (Detailed results will 
be shown in 0].) We therefore focus on the case with 
C = 2.47 x 10~ 3 . The typical value of P for our model is 
~ 10 34 erg/cm 3 (2.8M Q /M) 2 ; the initial maximum mag- 
netic field strength is then ~ 10 16 (2.8M Q /M) gauss at 
t = 0. (These scalings with M assume our adopted initial 
polytropic model, for which M/R = 0.22). This magnetic 
field is too strong to model a typical HMNS (but is sim- 
ilar in strength to 'magnetar' fields [lf|). However, the 
qualitative behavior obtained here still applies as long as 
the approximate scaling relation holds. 

Simulations in axisymmetry wereperformed indepen- 
dently using two GRMHD codes 7, 8], and the numerical 
results are qualitatively similar. As in many hydrody- 
namic simulations, we add a tenuous "atmosphere" to 
cover the computational grid outside the star. The at- 
mospheric rest-mass density is set to 10~ 7 p maX! o for the 
simulations shown here, where /3 maXj o is the maximum 
value of p at t = 0, which is 0.0668 in the adopted units. 

We perform simulations on a uniform grid with size 
(N,N) in cylindrical coordinates (w,z), which covers 
the region [0,4.5P] for each direction. For the HMNS 
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FIG. 1: Evolution of the central density, central lapse, and 
maximum values of \B X \ and \B y \ (the behavior of |B z | ma x 
is similar to the behavior of jB^jmax and is therefore not 
shown). IB 1 1 max and l-B^lmax are plotted in units of ^/pmax.o 
where p ma x,o is the maximum rest-mass density at t = 0. The 
solid (red), long-dashed (green), dashed (blue), and dotted 
(black) curves denote the results with iV=250, 300, 400, and 
500 respectively. The dot-dashed (cyan) line in the last panel 
represents the predicted linear growth of B y at early times. 



adopted here, R w 4.5M w 18.6 km(M/2.8M Q ). To 
check the convergence of our numerical results, we per- 
form simulations with four different grid resolution: 
A=250, 300, 400 and 500. We also checked that moving 
the outer boundary between 4P and 5R does not signif- 
icantly affect the results. 

Fig. ^ shows the evolution of the central density p c , 
central lapse a c , and the maximum values of |_B X |(= 
|P ro |) and \B V \{= n7|P v |) as functions of t/P c . Here 
P c w 39 M = 0.54(M/2.8M Q ) ms denotes the central 
rotation period at t = 0. The central density monotoni- 
cally increases with time up to the formation of a black 
hole. Evolutions with various grid resolutions demon- 
strate that the results begin to converge when N > 400. 
On the other hand, results are far from convergent for 
N < 300. For example, the maximum values of \B X \ are 
much smaller than those with higher resolutions, and the 
growth rate of \B X \ is underestimated. Hence, the effect 
of MRI, which is responsible for the growth of \B X \, is not 
computed accurately for low resolutions. This is because 
the wavelength of the fastest growing MRI mode is not 
well-resolved for low resolutions (see below) . 

Simulations with different treatments of the atmo- 
sphere are qualitatively the same when the atmospheric 
density is sufficiently small. However, the exact collapse 
time is somewhat sensitive to the details of numerical and 
atmospheric schemes. This is not surprising because, at 
late times, the star becomes marginally unstable, and 
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FIG. 2: The upper 4 panels show snapshots of the rest-mass density contours and velocity vectors on the meridional plane. The 
lower panels show the field lines (lines of constant As) for the poloidal magnetic field at the same times as the upper panels. 
The density contours are drawn for p/p ma x,o = I0~ a3l ~ 009 (i — 0-12). The field lines are drawn for A& = As >m in + (^^max — 
^40,min)j/2O (i = 1-19), where As,m&x and ^40, m i n are the maximum and minimum value of A^ respectively at the given time. 
The thick solid (red) curves denote the apparent horizon. 



the precise onset of collapse over the secular lifetime is 
sensitive to small differences in different schemes. 

For the chosen initial strength of the seed magnetic 
field, the early evolution is dominated by magnetic 
winding. When the seed field is weak, the induction 
equation shows that B y grows approximately linearly: 
B y (t;zj,z) « t-cuB l (0; zu, z)dift(0; zu, z). Indeed, the 
early growth rate agrees with the predicted one (cf. dot- 
dashed line in the last panel of Fig. . When the energy 
stored in the toroidal field becomes significant compared 
to the rotational energy, \B y \ grows more slowly and the 
degree of differential rotation is reduced. Eventually \B y \ 
reaches a maximum and starts to decrease. This is ex- 
pected to happen on the Alfven time scale , where 
the Alfven speed is va = \/b 2 /(ph + b 2 ). For the model 
considered here, the maximum value of va = 0.00734, 
and thus, the minimum value of tA — 15.8P C . We see 
that the maximum value of \B y \ starts decreasing when 
t > 2QP C , consistent with the Alfven time scale. 

MRI is evident at times t < 6P C as shown in Fig. ^ 
where the maximum value of \B X \ suddenly increases 
rapidly. MRI occurs wherever <9 ro S! < 19]. The 
wavelength for the fastest growing mode is Amri ~ 
2-kva/Q and the e- folding time of the growth is tmri = 
With our choice of the initial mag- 
netic field strength, Amri ~ -ft/10 and tmri ~ XP C . In 
Fig. ^ we see that MRI shows up prominently when N > 
400. Hence, we need to use a resolution A < Amri/10 



to study the effect of MRI accurately. We find that MRI 
first occurs in the outer layers of the star near the equa- 
torial plane. This can also be seen in Fig. [21 where we see 
that the magnetic field lines are distorted by t = 15P C . In 
Fig. ^ we see that the the central density begins to grow 
more slowly once \B X \ saturates. This may be caused 
by MRI-induced turbulence redistributing some of the 
angular momentum to slow down the contraction of the 
core. 

The combined effects of magnetic braking due to wind- 
ing and MRI eventually trigger gravitational collapse to 
a black hole at t « 66P C * 36(M/2.8M Q ) ms when an 
apparent horizon forms. The latest simulations [f| of 
BNS mergers show that for a sufficiently stiff EOS and 
typical observed BNS masses, HMNS formation is pos- 
sible. HMNS remnants are triaxial and strong emitters 
of gravitational waves in these simulations. The dissipa- 
tion time scale of angular momentum due to gravitational 
radiation is ~ 100 ms Therefore, HMNSs with an ini- 
tially large magnetic field (J5 > 10 16 G) will be subject 
to 'delayed' collapse due to MHD effects (magnetic brak- 
ing + MRI) rather than by the emission of gravitational 
waves. For seed magnetic fields which are much weaker 
than the cases studied here, gravitational radiation may 
be the trigger of collapse. However, it is possible that 
MRI may dominate the evolution even in this case, since 
the e-folding time of MRI is independent of the initial 
field strength. A more careful study of this scenario has 
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to be carried out in full 3D |20j- However, since any 
dissipative agent (viscosity, magnetic fields, gravitational 
radiation) serves to redistribute and/or carry off angular 
momentum, the final fate of an HMNS - collapse to a 
black hole, accompanied by a gravitational wave burst - 
is assured. 

Soon after the formation of the apparent horizon, the 
simulations become inaccurate because of grid stretch- 
ing. To follow the subsequent evolution, a simple excision 
technique for black hole spacetimes is employed 0, • 
The evolution of the irreducible mass of the black hole 
computed from the area of the apparent horizon Aah 
as Mi rr = ^/j4ah/167t, and the total rest mass outside 
the apparent horizon are shown in Fig. Soon after 
formation, the black hole grows rapidly, swallowing the 
surrounding matter. However, the accretion rate Mq 
gradually decreases and the black hole settles down to 
a quasiequilibrium state, i.e., the growth time scale be- 
comes much longer than the dynamical time scale. At 
the end of the simulation, M decreases to ~ 0.01M /P C . 
The estimated value of the black hole spin parameter is 
^lioic/^hoic ~ 0.8. The black hole angular momentum is 
computed from Jhoic = J — -/matter (r > Th), where the 
angular momentum of the matter outside the horizon 
•/matter (r > rjj) is computed by a volume integral (see 
e.g., Eq. (51) of 6J). The mass of the black hole Mhoic is 
crudely estimated from Mhoie ~ v^te + (^hoic/2Af irr ) 2 . 
The density contour curves and magnetic field lines at 
the end of the simulation are shown in the last column 
of panels of Fig. [3 

The value of Mq indicates that the accretion time scale 
is - 10-20P C w 5-10 ms(M/2.8M Q ). Also, we find that 
the specific internal thermal energy in the torus near the 
surface is substantial because of shock heating, indicat- 
ing that the torus can be a strong emitter of neutrinos. 



These facts suggest that the system formed after the 'de- 
layed' collapse of a magnetized HMNS (black hole + hot 
torus + collimated magnetic field) is a candidate for the 
central engine of short gamma-ray burs ts [23L 0] ■ This 
possibility is explored in more detail in |25j . 
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FIG. 3: Evolution of the irreducible mass and the total rest- 
mass outside the apparent horizon. 
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